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OVERVIEW 


This  report  summarizes  work  done  under  Contract  F046n-84-C-0032 
during  the  period  31  August  1984  through  1  March  1985.  The  principal 
goals  were  a  demonstration  of  state-space  procedure  for  modelling  and 
control  of  a  flexible  structure,  and,  more  generally,  increased 
understanding  of  control  problems  for  such  systems.  The  state-space 
approach  is  based  on  a  model  of  the  distributed  system.  That  is,  the 
model  includes  the  necessary  partial  differential  equations  wi thout 
modal  truncation.  The  basic  view  is  that  it  is  preferable  to  avoid 
introducing  approximations  until  they  are  required  (e.g.,  for  numerical 
calculations).  The  program  (formal  model,  state-space  model,  optimal- 
control  formulation,  approximation  procedure  and  numerical  calculation) 
is  carried  out,  in  detail,  for  a  simple  structural  system  including  a 
rotating  hub,  a  flexible  beam  and  a  tip  mass.  The  work  is  described  in 
Part  I  of  this  report.  A  final  section  in  this  Part  provides  some 
parallel  results  for  a  more  complex  structure  comprised  of  a  hub,  two 
flexible  support  beams  and  a  tip-body. 

The  second  part  of  the  report  is  concerned  with  some  parasitic 
effects  on  the  stability  of  a  distributed  system  with  feedback.  The 
approach  is  based  on  an  input-output  description  of  the  system; 
specifically,  a  transfer-function  approach  is  used.  In  order  to  keep 
the  calculation  burden  reasonably  small  the  structural  model  studied  is 
a  simple  cable  which  requires  only  second-order  spatial  oerivatives. 

The  baseline  system  employs  feedback  of  a  force  to  the  "free  end"  of  the 
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cable  at  a  magnitude  proportional  to  the  velocity  at  the  "free  end". 

The  gain  parameter  can  be  varied  to  produce  a  locus  of  roots;  albeit  one 
with  a  countably  infinite  number  of  branches. 

The  basic  system  is  examined  and  then  the  effects  of  three 
modifications  are  studied.  These  are:  time-delay  in  the  feed-back 
loop;  time-lag  in  the  feed-back  loop;  and,  viscous  damping  in  the  system 
forward-loop.  The  most  startling  result  is  that  with  time-delay  the 
undamped  system  will  be  unstable  for  any  non-zero  gain.  Viscous  damping 
apparently  provides  some  help  here;  with  such  damping  the  system  will  be 
stable  at  least  for  "small"  gain  values. 


•j 


INTRODUCTION 


This  retort  deals  with  modelling  and  control  of  the  structure  shown 
in  Figure  1.  It  is  assumed  that  the  structure  Is  pivoting  about  a  fixed 
pivot  at  point  0,  and  that  the  motion  is  in  a  plane.  The  structure 
consists  of  three  parts,  the  main  frame  (body  A),  the  beam  (body  B),  and 
the  end  mass  (body  C).  The  main  frame  is  assumed  to  be  a  rigid  body 
pivoting  about  point  0.  The  beam  is  assumed  to  be  flexible  and  rigidly 
attached  to  the  main  frame  as  a  cantilevered  beam.  Finally,  the  end 
mass  is  rigidly  attached  to  the  end  of  the  beam  so  that  it  moves  and 
rotates  witri  the  end  of  the  beam  and  its  attachment  point  is  assumed  to 
be  at  the  center  of  mass  of  the  end  mass.  A  tree  body  diagram  of  each 
mass  is  shown  in  Figure  2. 

This  report  may  be  logically  divided  into  five  parts.  In  the  next 
section  a  'formal'  model  for  the  system  Is  developed.  This  is  done  by  a 
direct  application  of  Newtonian  mechanics  to  each  of  the  structural 
components.  The  model  is  'formal'  in  the  sense  that  no  attempt  is  made 
to  show  that  the  resulting  system  of  (coupled  partial  and  ordinar  ) 
differential  equations  is  well-posed.  That  is,  we  do  not  prove  that  the 
system  has  a  unique  solution  for  a  certain  class  of  initial  data,  nor  do 
we  prove  that  the  solution  depends  continuously  on  the  initial  daca. 

A  second  section  is  concerned  with  development  of  an  abstract 
state-space  model.  Since  the  differential  equations  are  linear  (more 
appropriately  linearized),  the  abstract  model  is  in  the  form 
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y(t)  =  Ay{t)  +  Bu(t) 


where  >(t)  is  the  state  at  time  t,  u(t)  is  the  control  and  a  and  s  are 
linear  operators.  The  point  of  this  formulation  is  that  one  has  a 
theory  that  guarantees  wel 1-posedness  under  concrete  and  verifiable 
conditions  on  the  operators.  We  construct  our  model  so  that  the  system 
has  these  properties  and  hence  we  have  a  guarantee  that  the  differential 
system  makes  sense. 

A  third  section  is  devoted  to  linear-quadratic  regulator  theory 
for  the  abstract  system.  The  feedback  structure  of  the  optimal  control 
is  discussed  along  with  a  factorization  procedure  that  can  significantly 
reduce  the  computational  problem.  A  general  theory  of  tiumerical 
approximation  is  also  included  in  this  section.  A  fourth  section 
presents  details  for  developing  a  class  of  numerical  approximations. 

Some  numerical  results  are  included. 

Tne  final  section  describes  an  analysis  of  a  more  complex  structure 
including  a  rigid  hub,  two  flexible  side  beams  and  a  tip-body.  This 
analysis  parallels  that  of  the  simpler  case  but  provides  fewer  details 
and  no  numerical  results. 
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I.  EQUATIONS  OF  MOTION 


The  mathematical  model  of  the  structure  shown  in  Figure  1  is 
developed  in  this  section.  It  is  assumed  that  the  structure  is  pivoting 
about  a  fixed  pivot  at  point  0,  and  that  the  motion  is  in  a  plane.  The 
development  is  achieved  by  direct  application  of  Newton's  Laws  to 
individual  members  of  the  structure  with  the  resulting  equations  summed 
to  determine  the  overall  motion.  The  development  is  rigorous  with 
assumptions  stated  as  needed  or  desired. 

The  structure  consists  of  three  parts,  the  main  frame,  or  mass  A, 
the  beam,  or  mass  B,  and  the  end  mass,  or  mass  C.  The  main  frame  is 
assumed  to  be  a  rigid  body  pivoting  about  point  0.  The  beam  is  assumed 
to  be  flexible  and  rigidly  attached  to  the  main  frame  as  a  cantilevered 
beam.  The  end  mass  is  rigidly  attached  to  the  end  of  the  beam  so  that 
it  moves  and  rotates  with  the  end  of  the  beam.  Finally,  the  attached 
point  is  assumed  to  be  at  the  center  of  mass  of  the  end  mass.  A  free 
body  diagram  of  eacn  mass  is  shown  in  Figure  2. 

Beam  Equations 

The  coordinate  system  to  be  used  is  fixed  in  the  main  frame  and 
rotates  with  it.  The  origin  is  at  point  0  with  the  x  axis  pointing 
along  the  undeflected  beam,  the  z  axis,  the  axis  of  rotation,  and  the  y 
axis  completing  the  right  hand  set.  The  position  vector  to  some  point 
on  the  beam  is  given  by 
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(1. 


r(t,x)  =  r(x)  +  c(t,x) 


where 

r(x)  =  position  to  undeflected  location  of  mass  element  dm 

Uit.x)  =  position  of  mass  element  dm  with  respect  to  its 
undeflected  position. 

It  IS  assumed  here  that  the  deflections  U  are  small.  The  governing 
equations  are  given  by 


r(t,x)  »  r(x}  +  u(t,x) 

^rel^t.^)  .  -'•'u(t,x]  ^  ujt.x) 

ot 

where 

n  I't.x)  =  velocity  relative  to  the  rotating  coordinate 
system 

r^  l(t,x}  '  acceleration  relative  to  the  rotating  coordinate 
system 

and  by  Newton's  Laws  which  are  [6]; 

F  -  '  d  dm 
'  m 

~  .  r  X  a  cm 
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where 


F  -  total  external  force 

Mq  =  total  external  moment  acting  about  point  0. 

The  general  expression  for  the  acceleration  at  any  point  is  giver  by  'c 


.-X  (.X  r)  ^x  r  +  2.,  X  r 


he  can  note  that  for  the  conditions  of  this  problem  7^-0  since  it  is 
fixed  pivot. 

We  can  apply  the  moment  equation  (1.3)  to  the  beam.  Substituting 
the  expressions  for?  and  its  derivatives  (eq.  1.2)  along  with  the 
expression  for  the  acceleration  (eq.  1.4)  into  (eq.  1.3)  yields  the 
general  result: 


=  'n  '■  r  [“  X  (^-  X  r) 


/]  +  u  X  [>■•  X  ('X.  X  r)]+  r  x 


X  (j  X  u)j  +  u  X  X  (ju-  X  u)]  f  r  X  (  ...  X  r)  +  u  X 


^  _  _ 

(..  X  r)  +  u  X  (ju  X  u)  +  r  X  x  u)  -*■  2r  .<  (  .  x  u ^  ■ 

»  2u  X  (..  X  u^)  +  r  X  Uj.^  u  X  •  dm  (l-S/ 

—  f  P  ' 

wiiere  Mg'''  is  the  total  applied  mcment  on  the  beam  about  point  0. 
Equation  (1.5)  reduces  considerably  if  we  take  advantage  of  some 
.additional  assumptions.  If  .  is  assumed  small  so  that  products  such  as 
■  .-u  and  0^  course  j'  are  negligible,  all  but  two  terms  nr-e  dropped 


from  eq.  (1.5).  If  we  further  assume  the  deflections  occur  only 
perpendicular  to  the  beam  (no  stretching  of  the  beam,  consistent  with 
small  angular  rates  and  deflections),  then  u  »  u j ,  and  eq.  (1.5)  oecomes 


(lQga,(t) 


m  tt 


(t,x)  dm)  k 


(1.6) 


where 


‘OB 


=  X  (T  X  r*)  dm 

=  moment  of  inertia  cf  the  undeflected  beam  about  point  0. 


Equation  (1.6)  is  the  beam  moment  equation. 

The  elastic  equation  for  the  beam  is  obtained  from  basic  prin¬ 
ciples.  If  the  beam  is  assumed  to  have  uniform  properties  along  its 
length,  and  that  no  external  forces  act  on  the  beam  except  those  at  the 
boundaries,  then  the  equation  of  motion  of  an  element  of  the  beam  is 
given  by  [5) 


ait,x: 


=  .  il  -u(t.x) 


J  =  -  — 


El  u  (t.x)  j 
xxxx  ^ 


JX 


wherG 


1,7; 


El  =  the  stiffness  properties  of  the  beam  (constant) 

-  density  per  unit  length  of  the  beam  (constant) 

T(t,x)  »  the  acceleration  of  the  beam  element. 

The  acceleration  is  given  by  eq.  (1.4)  which  can  be  substituted  into  eq. 
(1.7).  Under  the  previous  assumptions,  eq.  (1.7)  reduces  to 


u^^(t,x)  =  -n  u^^^^(t,X;  -x...(t,- 


(1.8) 


G 


Main  Frane 


The  equations  for  the  main  frame  are  sin’’'  the  rigid  body 
equations  of  motion.  The  moment  equation  is 

^(t)  =  kt)  k  (1.9 

where 

=  moment  of  inertia  of  main  frame  about  point  0. 

End  Hass 

Both  the  force  and  moment  equations  are  needed  for  the  end  mass 
since  it  is  translating  as  well  as  rotating.  The  force  equation  is 
given  by 


^  \  “  mc7(t.L) 


(1.10) 


where 

f|  =  force  on  mass  due  to  beam. 

Under  the  previous  assumptions,  eq.  (1.10)  becomes 

^  (1,11) 

where  n(t'  =  u^(t,L),  the  relative  velocity  of  the  end  cf  the  beam 
(see  eq.  1.2). 


7 


For  planar  motion,  only  the  z  component  of  the  moment  1s  of 


interest  and  leads  to 


Kat) 


(I.IZ) 


where 

I  =  moment  of  inertia  of  the  end  mass  about  its  center  of 
mass 

<+>(t)  ■  the  angular  position  of  the  end  mass  with  respect  to  an 
inertial  reference. 

The  angular  position  of  the  end  mass  is  that  due  to  the  rigid  body 
motion  of  the  beam,  (5),  plus  the  angular  motion  due  to  the  deflection 
of  the  beam  which  is  equal  to  the  slope  of  the  beam, 

d{t)  =  e  (t)  + 

lit)  =  +  ;(t)  (1.13) 

V(t)  =  :(t)  +  ;(t) 


where 


^(t;  =  u^^(t,L,],  the  angular  rate  of  the  beam  end  due  to 
def 1 ection . 


Consequently  the  moment  equation  can  be  written  as 


M 


C 


(t;: 


•  *‘+. 
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System  Moment  Equation 

The  moment  equation  for  the  complete  system  can  be  obtained  by 
combining  the  moment  equations  for  each  individual  mass.  In  order  to 
carry  out  this  step  it  is  necessary  to  write  the  expressions  for  the 
applied  moments  on  the  main  frame  and  on  the  beam.  These  are 
respectively  (see  Figure  2) 


^0^^'  ^  ■  ^08 


/ 1 


1.  is; 


and 


(B) 


"08 


M..  *  LFi 
c  Ly 


(1.16) 


By  adding  all  the  moment  equations  (1.6,  1.9  and  1.14)  together  and 
utilizing  eqs.  (1.11,  1.15  and  1.16)  the  system  moment  equation  can  be 
obtained.  The  result  Is 


^OB  *  ^  m^L^)  •u)(t)  +  *  ni^Li(t) 

L 

+  /  ^  ,X/cdx  (1.17) 

0 

where  the  relation  dm^odx  has  been  used. 

Boundary  Condi tions 

Associated  with  the  structural  equation  (1.8)  are  the  following 


boundary  conditions;  (see  Figure  2) 


At  X  ■  0 


u(t,0)  >  0 


Uj^(t,0)  »  0, 


"oe  ■  -E' 


cantilever  moment  to  mainframe 


(1.18) 

(1.19) 

(1.20) 


At  X  =  L 


u^(t.L)  »  n(t) 


(1.21) 


u^^(t,L)  »  ^ (t) 


(1.22) 


El 


(1.23) 


(1.24) 


Final  Working  Equations 

For  the  analysis  which  follows  in  the  next  section,  It  is  useful  to 
rearrange  the  equations  developed  previously  into  a  final  working  set. 

We  have  the  structural  equation  (1.8) 


u^^(t.x)  ^  xAt)  =  -El  u,,,^(t.x). 


(1.25) 


’0 


From  the  main  frame  moment  equation  (1.9),  eq.  (1.15)  and  the  boundary 
condition  (1.20)  we  obtain 

I^<L(t)* -El  u^j^(t,0)  +  M^(t)  .  (1.26) 

The  end  mass  force  and  moment  equations  (1.11  and  1.14  respectively)  in 
conjunction  with  the  boundary  conditions  (1.23)  and  1.24)  become 

I<.[w(t)  ^  ;  (t)]  »  El  u^j((t.L)  (1.27) 

and 

rncCL<L(t)  +n(t)]  «  El  u^^^it.L)  (1.28) 

Equations  (1.25-1.28)  plus  the  remaining  boundary  conditions,  eqs. 
(1.18,  1.19,  1.21  and  1,22)  are  the  starting  point  for  the  analysis  in 
the  next  sections. 
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II.  STATE-SPACE  FORMULATION 

The  purpose  of  this  section  is  to  construct  a  state-space  model  for 
the  system  dynamics.  That  is,  we  are  seeking  an  appropriate  state-space 
Y  so  that  the  (linearized)  equations  can  be  written  as 

y(t)  «  Ay(t)  +  flf(t)  (2.i; 

where  f  is  the  control  and  a  and  s  are  appropriate  operators.  In 
addition  we  shall  show  that  the  operators  a  and  b  are  sufficiently 
"nice"  that  the  system  (2.1)  is  well-posed.  This  means  that  for 
appropriate  initial  data  [y(0)j  and  control  input  f,  the  system  (2.11) 
has  a  unique  solution  and  that  solution  depends  continuously  on  the 
Initial  data. 

The  formal  modelling  procedure  of  the  previous  section  resulted  in 
the  following  description: 


U„(t.x)  .  .  -a  U„„(t,x) 

D 

(2.2) 

I^^t)  «  El  Uj^x^t.O)  +  M^(t) 

(2.3) 

I^['j(t)  -t-  C(t)]  •  -El  u^^(t,L) 

(2.4) 

mj.[L‘^(t)  +  n{t)]  •  El  u^j^^(t,L) 

(2.5) 

12 


In  addition  to  these  differential  equations  there  are  important 
boundary  conditions.  Since  the  beam  remains  Joined  to  the  hub  one  has 


u(t.O)  »  u^(t,0)  ■  0 


The  cantilever  nature  of  the  connection  also  requires  that 


u^(t,0)  *  «  0  , 


while  integrity  of  the  upper  joint  requires 


u^(t,L)  »  n  (t)  , 


and 


Uj^^(t,L)  =  C(t)  . 

As  a  first  step  in  the  state-space  formulation  observe  that  the 
partial  differential  equation  (2.2)  can  be  re-cast  as  a  first-order 
system. 

5[U  (t,x)]  *  Jj.^[u  (t.x)  +  X'.(t)] 

It  Ix^  ’■ 

(2.6) 


J_[u  (t,x)  f  x..^(t)]  =  -El  (t.x); 

St  Sx2 
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This  suggests  that  the  functions  (of  x'  + 

x*'jj{t)]  are  worthy  of  consideration  as  state  components.  With  this 
insight  define  a  "state"  2(t)  by 

aj(t) 
n(t) 

2(t)  -  4(t) 

"xx"-’'' 

U^(t,x)  +  xw(  t) 

Note  that  z(t)  has  five  coordinates: 

■oi  angular  velocity  of  the  hub 

n  relative  lineal  velocity  of  the  tip 

?  relative  angular  velocity  of  the  tip 

u^^  curvature  of  the  beam 

u*+  X  velocity  of  the  beam 

t  U) 

The  quantities  ''  and  C  are  each  scalar  while  3nd  [u^  +  x 

are  functions  of  the  spatial  variable  x.  Thus,  the  state-space  can  be 
at  least  formally  identified  as 

Z  =  R  X  R  X  R  X  X 

where  R  denotes  the  real  line  and  ^9  is  the  usual  Lebesque  space  of 
real-valued,  square-integrable  'functions’  defined  c  1  the  interval  [0, 
(see  [9]). 
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It  1s  possible  to  re-write  the  equations  (2.2)  -  (2.5)  in  terms  of 
the  five  components  of  Z  and  the  control  M^. 

ij(t)  *  El  z^{t)  *  M^(t) 

[L  z^(t)  +  2^(t)]  =  £I  ^  oDZ^(t) 

(2.7) 

C2^(t)  +  Z2(t)]  =  El  24(t) 

Z^(t)  *  25(t) 

Z5(t)  -  tl  0^  2,(t) 

u 

In  these  equations  d  denotes  differentiation  with  respect  to  the  spatial 
variable  x  and  denotes  evaluation  at  x  »  b.  Boundary  conditions 
require  that 

ZgCt)  =  u^(t,0)  =  0  , 

^5^^^  *  u^(t,L)  +  L4t) 

=  r,{t)  +  Laj(t) 

»  z^Ct)  +  L  z^(t)  , 

ffQoDZjCt)  =  Uj.^it,0)  =  0  , 
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and 


oDZ^(t)  *  U^^(t,L)  +  ju(t)  .  L 

=  z^it)  +  z^(t)  .  L. 
The  system  (2.7)  is  in  the  form 

f2(t)  =  CTZit]  +  //f(t)  . 


In  principle,  this  can  be  brought  into  the  normal  form  (2.1)  by 
'inverting'  r.  It  is,  perhaps,  somewhat  easier  to  introduce  a 
coordinate  transformation  that  'uncouples'  tne  left  hand  side. 

Let  y.j  *2^  (i  *  1,  4  and  6)  while  y^  *  LZj  +  and  y,  = 

2j  +  verified  that  the  system  (2.7)  can  be  written  in  the 

form  (2.1)  with 


.<1  = 


and 


0  0 

C  0 

0  0 

C  0 

0  0 


0 

0 


n 

Ia 

EI 


"O 


D 

0 


0  _  EI  0 


.e; 
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The  boundary  conditions  are  readily  translated  to 

yg(t)  ~ 

o^ys^^)  "  y3(t) 


These  conditions  are  incorporated  by  restricting  the  domain  of  the 
operator  a  to 

o(-'.)  =  {(yp  y2.  y3.  y^)  c  x  x  j 

y4.  ys  £  y5(0)  =  O,  ygd)  =  :2.1 

y5(0)  --=  0,  y5(L)  =  y3'  . 

2  2 

W  ’  is  the  u-ual  space  of  real-valued  functions  with  the  function  drc 
its  derivative  ini:^[l  ].  As  a  distance  measure  on  the  space  Y  we 
introduce  the  inner-product 
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<>.W>  =  \  ^  V2'"2  "  i  ^3"^ 


+  El  yAx)  w.(x)  dx  +  /  *- 
0  ^  0 


>5' 


dx 


C.i: 


It  can  be  verified  that  ly';  is  a  norm  and  that  Y  is  q  Hilbert 

space.  The  inner-product  ^y,y>  is  the  mechanical  energy  in  the  physical 
system  at  state  y. 

At  this  point  we  have  constructed  a  system 


y(t)  =  Ay(t)  +  flf(t)  (2.1) 

with  state-space  Y  »  x  x  inner  product'^*  [given  by  (2.111[ 
and  operators  ^  and  ®  [given  by  (2.8),  (2.9)  and  (2.10)].  The  final 
part  of  this  section  Is  concerned  with  establishing  that  this  system 
makes  sense.  Formally  one  can  write  down  the  solution  to  (2.1)  in  terms 
of  the  variation  of  constants  formula 

y(i.)  =  e^S(O)  f(s)ds 

Since  3  is  bounded,  the  solution  y  will  have  the  desired  orooerties 
(existence,  uniqueness  and  continuous  dependence  on  initial  data;  if, 
and  only  if,  e"^^  is  a  Cq -semigroup .  Hence,  we  niust  show  that  an  a 
operator  "generates"  such  a  semigroup.  This  is  the  central  matter  in 
abstract  formulaf  ons  such  as  (2.1). 

Over  .  5  u  thirty-five  ^ears  there  has  been  significant  progress 

in  charart  -ii'  :u  ‘hose  operators  that  are  generators  of  semigroups 
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[4,7,12].  We  shall  make  use  of  a  special  form  of  the  Lumer-Ph i 1 1 i ps 
Theorem  [7]. 

Theorem  (Lun.cr-Phi  1 1  ips) ;  Let  ^  be  a  densely  defined, 
closed  linear  operate-.  If  both  «  and  .a  are  dissipative, 
then  A  is  the  generator  of  a  Co  semigroup  of  contractions 
on  y. 

To  employ  this  theorem  we  must  demonstrate  that  the  «  defined  by 
(2.8,  2.10)  satisfies  the  hypothesis.  The  fact  that  is  closed  and 
densely  defined  follows  from  standard  results  in  theory  [9], 


After  elementary  calculations  it  1s  seer  that  each  term  has  a  common 
factor  (|.EI).  The  last  two  terms  are  each  integrated  by  parts  to 
produce 

-  y^y4(0)  +  y^y^d)  -  y3y4(L) 

+  ty4(L)y5(L)  -  y4(0)y^(0)  -  ^  y^(x)  y^(x)  dx; 

-  (y5(L)  y^d)  -  y5(0)  y'^{Q)  -  /'-  y^(x)  y^(x)  dx}  . 

The  integral  terms  cancel  and  the  remaining  terms  can  be  gathered  to 
yield 

-  y4(0)[y^  -  y^(0)]  +  y^d)  [y^  -  ygd)] 

+  y4(L)  Cy^d)  -  y3]  *  y^(0)  [y^{0)] 

It  is  now  easily  verified  that  for  y  in  the  domain  of  [see  (2.11)1  each 
of  the  bracketed  terms  is  zero  so  that 

<  y,  ^y  '  =  0  . 


Note  that  we  have,  in  fact,  shown  that  the  operator  is  conser¬ 
vative.  Since  our  norm  is  related  to  the  energy  we  have  simply  verified 
that  the  energy  is  constant  {for  the  uncontrol  lea  system).  Tc  complete 
the  hypothesis  of  the  theorem  we  most  verify  that  the  adjcint  operate’-, 4* 
is  also  dissipative.  In  fact  we  siiall  show  more. 
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;.emma  1.2.  The  operator/;  is  skew-adjoint. 


Proof:  Recall  that  the  operators  is  defined  on  its  domain 

*  *  * 

d{a  )  =  {w/ there  exists  w  e  2  suc^'  that  <Ay,"/''  -  <y,w  ■  =  0 

for  all  ye  d(a)) 

by 

»  * 

A  W 

We  use  this  characterization  to  compute  «*  . 

Let  y  *■  col  (yj  ,y2»y3 >y4( • )  •  ))e  assume  w  » 

col(Wj,W2,W2,w^(.),W5(.))and  w*  ■=  col  (wj*.  ,  w^*,  w*(.),  w*(,))  satisfy 

♦ 

0  =  <«y ,w?  -  ^y.w  > 

Using  the  definition  of/;  and  the  inner  product  <  ,>  we  obtain  the 
identi ty 

0  =  [y4(0)(EIw,)  *  y4  (L)  (EIw^)  -  y4(L)  (Elwj) 

+  El  -'^>c(x)w.(x)dv  -  EI/^  y''(x}wi.(x)dx  ] 

0  ^  ^  0  ^ 

-  [yi(lAWi) 

El  A  y  ,fx;w*(x  )dv  +  yjx)w*(x)dv  ]  . 
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I  I 

Collecting  terms  and  applying  the  conditions  that  y^  »  >'5^0',  y-.  =  y^  (L), 
y^  =  >'5(1')  and  yg(0)  »  0,  this  equation  becomes 

0  *  y^yO)  +  y^CL)  (EIw^)  +  y4(L)  (-EIW3) 

^  ^5(0)  ^■‘OA'^P  *  ■*•  y5(L )  (-Ij.w*) 

+  [  y^(^)  (£Iw^(x;)  +  yjfx}  (-ow*[x))]  dx 

C  y^i)^}  (-EIw^(x))  ■►y^fx)  (-EIw^(x) dx, 

which  must  hold  for  all  functions  y^j'x),  ye(x)  in  the  Sobolev  space 

''^(0,L)  (see  [1]).  In  particular,  this  equation  hclds  for  y.(x)  *  0  and 

y^(x}  satisfying  y.(0)  »  ylCo)  =  y^d)  =  ydL)  =  0,  i  ,j  *4,5  which 
O  J  J  J  J 

implies  that 

y'^fx)  (EIw^(x))  +  yjfx)  (-0V^(x;)]  dx  =  0 
and 

q'"  y4(x)  (-Ei'w5(x))  +  y^(x)  (-EIw*(x))]  dx  *  0. 

The  fundamental  Lemma  of  the  Calculus  of  "ari  •itio'^s  innlies  that 

2 

w^i'.;  and  Wg(.)  belong  to  .v  (0,L)  and 


I 

i 

i 

k 


and 

|I  w^'fx)  *  w*(x). 

Substituting  for  w5(y)  and  w|(x)  in  the  Integrals  and  Integrating  b> 
parts ,  the  equation 

0  =  <y,w>  -  <y,w> 

becomes 

0  =  y4(0)  [Elwj  .  £lw^(0)]  ♦  y5(0)(:£Jw;^(0): 

*  y^(0)  [EIw^(0)J  +  y^(0)  [-£lw^(0)  -  I^h*] 

*  ^'4^'-^  ^1^(1-)]  +  y5(L)  [.EIw^jfL)-  wf 

^  >'4^1-)  [EIw^  -  £iWj(L)]  f  y^fL)  [£Iw^fL,'  - 

which  must  hold  for  all  values  of  y^(0).  y,(0) .  y.d; .  y,(L).  y^^'O), 
ysiL).  yjd)  and  for  y^iO]  =  o.  Thus,  it  follows  that 

WgfO)  =  w^,  w^(L)  =  w^,  w;(L)  »  W3,  w,(0)  »  0 
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and 


w*  =  -El  w^(0),  ^2  »  -  EJ[  w‘^(L),  w*  =  El  w^(L)  . 

'"c 

# 

Therefore,  if  w  e  d(a  ) ,  then  w  *  col (w^ ,  w^,  . ) ,  Wg( . )  is  such 

2 

that  w^(x),  w^(x)  e  h  (O.L)  with  Wj  ^  W2(0),  =  w^d)  ,  ,  0 

=  w^(0) ,  and 

/I w  -  col  (-^  w.(0) ,  •  JJ  .  'WcCx)  ,  £I  w.(x ; ;  * 

'a  “c  'c  '' 

Conversely,  it  1s  easy  to  show  that  if  w  is  as  above  then  wco(a  )  and 

*  0 

w  *  -''w.  This  completes  the  proof  that  o  {a  )  =  d{  a)  and 

a\  =  -  a\h  In  particular,  -a  Is  skew-adjoint. 


Combining  Lemmas  2.1  and  2.2  we  have  established  the  following. 


Theorem  2.1:  The  operator  defined  by  equations  (2.8)  and 

(2.10)  generates  a  Co  semigroup  of  contractions  on  the  state- 

3 

space  Y  =  R  X  ^2  * 

We  have  verified  that  the  state-space  model  assembled  from  the 
formal  system  does  indeed  make  sense.  Before  closing  this  section  we 
note  that  the  angular  coordinate  0  may  be  appropriate  in  the  control 
problem.  That  is,  we  want  to  add  a  "zeroth"  state,  with  equation  I(t;  = 
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To  avoid  duplicating  much  of  the  previous  work,  we  shall  do  this 
augmenting  the  state  and  "perturbing"  the  A  operators. 

Let  2  ”  R  X  Y  and  consider  the  dynamical  system 

z(t)  =  A  2(t)  +  flf(t) 

wi  th 

A  (a.y)  = 

(2.12) 

B 

and 

d{a)  =  {  (a.y)  e  R  X  Y|  cx  c  R  and  y  £  ri{A)  >  . 

The  inner-product  is  extended  in  an  obvious  way 
(ii.w)  z  *  <y.w>y  • 

It  is  clear  that  one  can  write  as  a  sum 


where 


-^0 


0 

Ay 


and 


•^2  »y )  “ 


-i5" 
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It  can  be  easily  seen  that generates  a  contraction  semigroup  on  2. 
The  proof  is  exactly  the  same  as  that  used  to  demonstrate  that  ^ 
generates  a  contraction  semigroup  on  Y.  Also  it  is  clear  that is  a 
bounded  linear  operator.  A  standard  result  [7,  p.  7G]  can  be  used  to 
conclude  that  a  generates  a  Co  semigroup  on  Z.  It  should  be  noted  that 
A  is  not  conservative  and  hence  ||2(t)j|  will  not  be  constant  in  the 
unforced  case,  b  is  obviously  bounded. 


III.  THh  CONTROL  PROBLEM 


In  this  section  we  formulate  a  general  quadratic  optimal  control 
problem  for  distributed  parameter  systems  and  then  review  those  results 
that  are  important  for  developing  approximation  and  numerical  schemes 
for  computing  feedback  gains.  Let  Z,  U  and  A  be  real  Hilbert  spaces. 

We  use  the  symbol  <,>  to  denote  the  inner  product  in  each  space.  The 
space  Z  is  the  state  space,  U  is  the  space  of  control  Inputs  and  A  is 
the  space  of  outputs  (observations),  .‘furthermore,  we  assume  that  a  with 
domain  d{a)^  I  is  the  generator  of  a  strongly  continuous  semigroup  e'^^ 
on  Z  and  b;  U  -►  Z  and  v.  Z  -►  A  are  continous  linear  operators.  We 
assume  that  r:  II  -*■  U  is  continuous,  linear,  self-adjoint  (i.e.  r  *  r*) 
and  1 1.‘?|  I  ^  m  >  0.  Moreover,  we  define  Q  »  v*v 

The  distributed  parameter  control  system  is  governed  by 

z(t)  »  /2 (t)  +  su(t)  (3.1) 

with  output 

y(t)  -  vzU)  .  (3.2) 

Given  z^,  the  (mild)  solution  to  (3.1)  satisfying  z  (o)  is  given  by 

the  variation  of  parameters  formula 

z(t)  =  e-’^z.  +  e^U-s)  ^ 

^  0 
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We  consider  the  problem  of  finding  u(-)  c  12(0, T;U)  which  minimizes 

J  (u)  «  {<  Q  z  (s),  z(s)^+  *" /?  u  (s),  u(s)  ->)  ds  ,  (3.4) 

0 

where  z(t)  is  given  by  (3.3)  for  01  tl  T.  This  problem  is  the  most 
direct  extension  of  the  standard  finite-dimensional  linear  quadratic 
control  problem. 

It  can  be  shown  (see  Gibson  [3])  that  the  optimal  control  u*’(t)  is 
state  feedback.  In  particular, 

u*(t)  -  -  5*n(t)z(t)  ,  '3.5) 

where  n(t)  :  Z  Z  is  the  solution  to  the  operator  Riccati  integral 
equation 

.T(t)  -  e^'  iQ  -  n{r,)  fl*n(n)]  dn  .  (3.6) 

t 

It  is  important  to  note  that  in  general  (3.6)  cannot  be  differentiated 
to  produce  a  Riccati  differential  equation.  This  lack  of  "smoothness" 
is  one  reason  that  numerical  methods  for  computing  approximate  (sub- 
optimal)  gains  have  to  be  carefully  designed. 

Another  feature  that  is  apparent  from  (3.6)  is  that  the  adjoint 
operator  a*  plays  an  important  role  In  determination  of  :^(t).  For 
infinite-dimensional  problems  such  as  structural  control,  the  operator 
is  a  differential  operator.  It  is  not  a  continuous  operator.  The  same 


is  true  of  a*  .  approximation  scheme  used  for  computing  sub-optimal 
gains  must  be  based  on  a  numerical  algorithm  that  produces  good 
approximations  to and  a  .  This  turns  out  to  be  a  crucial  point  that 
is  often  overlooked  in  finite  element  and  modal  control  approaches  to 
these  problems. 

As  we  shall  see  below,  computing  numerical  solutions  to  the  Riccati 
equation  (3.6)  involves  solving  non-linear  mixed  integral  and  partial 
differential  equations.  This  is  a  highly  complex  problem.  There  are 
other  methods  for  computing  the  gain  operator 

K(t)  »  B*n(t) 

directly.  One  such  method,  the  Chandrasekhar  algorithm,  is  ideally 
suited  for  these  infinite-dimensional  problems.  In  particular,  it  can 
be  shown  (see  [8])  that  the  gain  operator  can  be  found  by  solving  the 
coupled  Chandrasekhar  integral  equations 

K(t}  =  s't'  (nlKnld-.  (3.7) 

t 

:(t)  =  :,(-)  b  K{r}  d-  ,  (3.8) 

t 

where  A(t):  Z  -  U  and  r,{t):  Z  -  A  are  bounded  strongly  continuous 
operators . 

In  many  applications  equations  (3.7)  -  (3.8;  are  less  difficult  to 
solve  numerically  than  equation  (3.6).  This  is  particularly  true  for 
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problems  with  a  fixed  finite  number  of  sensors  and  actuators.  Such 
problems  occur  in  realistic  formulations  of  structural  control  problems, 
and  result  in  finite  dimensional  control  and  observation  spaces,  i.e. 

U  and  A  “P^. 

In  order  to  illustrate  the  potential  advantages  of  the 
Chandrasekhar  algorithm,  consider  the  control  problem  for  the  structure 
described  In  Sections  1  and  2.  The  state  space  is  Z  =  P*^  x  (0,L) 

X  1’^’®  Riccati  operator  n(t}:  Z  ^  Z  has  the  form 


where 


00 


(t)  Is  a  4  X  4  matrix  and  ii^.(t)  have  the  representations 


'01 


(t) 


(.) 


:2(.) 


/ 


Pl(t.T) 

nJlt.T) 
r  (t,T) 
s}(t,T) 


P2 


(t,: 

2 


ro(t,0 
S2  ( !•  >^  ) 


0  (  T  ) 

n 


’"'1  ■ 

(x)  = 

Pl(t,x) 

qi(t,x)  r^(t.x)  s^(t,x) 

"1 

'^2 

"3 

q2(t,x)  r2(t,x)  S2(t,x) 

-  "'4  - 
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and 


where  ,  q . ,  r . ,  i ,  j  «1,2  are  functions.  The  Riccati 

equation  (3.6)  becomes  a  coupled  Integral -partial  differential  equation 
that  must  be  solved  to  determine  the  16  entries  IniiQ^it)  and  the  12 
functions  p.,  q . ,  r . ,  s.  anda^j  for  1.j-l,2,  0<  t^  T,  0^  x,_^  L. 

On  the  other  hand  since  the  operator  b  >  col  l0,l/I^,0,0,0,0) 
has  rank  one  the  gain  operator  ''‘(t)  ■  «  ^fi*(t)  has  rank  one  and  /((i,  has 
the  form 

^  i.  L 

/<(t)z(t)  =  :  k.(t)2.(t)  +  /  K.(t,s)z,.(t.s)ds  +  /  K,(t,s)2.(t,s)ds 

1=1  ’  ’  0  ^  ^  0  ®  ® 

or  equivalently, 


^(t)z(t)  =  kj(t)e(t}  +  K,(t)  -^(t)  +  K^(t)  -  (t)  +  K^{t) 


+  /'  kj(t,s)u^^(t.s)ds  Kg(t,s}[u^.(t,s;  ;os  (3.9) 


In  order  to  determine  the  gain  operator  it  is  sufficient  to  determine 
the  6  functions  ,  i  =  1,6.  The  functions  Kc(t,s)  and  Kg(t,s;  are 
called  the  functional  gains.  These  functions  are  linear  combinations  of 
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P^.  >■.,  s . ,  1  =  1,2  are  not  needed 


q^,  (t,s).  Note  that  the  functions  -:'  ^.,, 
to  define  the  gain  operator. 

p 

If  the  observation  operator’*'  has  finite  rank,  say  :  Z  ^  . 

then  -  (t)  has  the  form 

4  ^ 

I{t)z=  I  A.(t)2,-  Ag(t,s)Zc(s)ds  + /'-  A,'t,s;z,(s}ds  , 

where  the  functions  A.(t)  take  values  in  i  =  1,2, 3, 4  and'A.Jt.s)  : 
for  i  »  5,6.  In  this  case  the  Chandrasekhar  equations  involve  coupled 
integral -partial  differential  equations  to  be  solved  for  the  6  functions 
K.j(t,s),  0  t^  T,  0^  s^  L  and  the  6 'p  functions  that  are  the 
entries  in  the  valued  vector  functions  A,j(t),  A.(t,S/  for  01  t^l  T, 
01  si  L.  The  main  point  to  be  made  here  is  that  the  spatial  dimension 
of  the  equations  (3.7)  -  (3.8)  is  one,  compared  to  equation  (3.6)  which 
has  spatial  dimension  two.  This  reduction  can  be  substantial  when  one 
attempts  numerical  solutions  cf  the  equations. 

Perhaps  the  simplest  way  to  see  the  real  advantage  of  the 
Chandrasekhar  algorithm  is  to  apply  it  directly  to  the  fi’ite  dimen¬ 
sional  approximating  system.  An  example  cf  this  type  will  be  given  in 
Section  4. 

In  order  to  obtain  approximate  solutions  to  the  Riccati  equation 

(3.6)  or  the  Chandrasekhar  eouations  (3.7)  -  (3.8)  it  is  necessary  tc 

■*’  ■* 

approximate  the  operators  ;i  ,  a  ,s.p,  v  and  v*  .  ’’he  most  difficult 
aspect  of  this  problem  is  the  development  of  convergent  apprcximaf or 
schemes  for  .a  and  «*.  We  shall  need  the  following  definitions. 
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Cefinition  1.  An  approximating  sequence  for  the  concrol  problem  defineo 
by  equations  (3.3)  -  (3.4)  is  a  sequence  {«  ,  ®  )  such  that  a'  , 

andC^ 


v''*  are  continuous  linear  operators  satisfying 


j  R^\  I  il  '''  1  0>  is  non-negative  definite,  •*  fiu  and  ^^'u  ^^u  for 

N  N 

dll  u  -  U,  ^  2  w  z  for  all  2  ^  Z  and  the  operators  «  generate  strongly 

continuous  semigroups  e  ^  satisfying 


N 

At. 


e^^2 


(3.1 


for  all  2  s  z,  uniformly  for  te  [0,T] 

An  approximating  sequence  is  said  to  be  a  strong  approximating  sequence 

r^^'t  r  Ni*  .  , . 

If  the  semigroup  e^'  ■'  generated  by  [a  j  satisfies 

e  -  e^*'w  (3.1Z 

for  all  w  Z  uniformly  for  t  :  [O.T^ 


Ine  construction  cf  approximating  sequences  for  the  control  problem 

(3.3;  -  (3.4;  is  a  problem  in  numerical  analysis  and  approximation 

theory.  The  basic  idea  is  to  approximate  the  (differential)  cperator 

N 

by  a  finite-dimensional  operator  a  (i.e.  using  finite  elements,  finite 

N 

differences,  modal  truncation,  etc.)  and  then  showing  that  a  converging 

to  A  implies  (3.11)  £nd  (3.12).  This  is  a  nontrivial  problem  in 

N 

functional  analysis.  In  fact,  the  basic  question  (when  does  az 

imply  (3.11)  and  (3.12))  is  not  yet  fully  understood.  A  partial  answer 


to  the  convergence  question  is  provided  by  the  now  famous  Trotter-Kato 
Theorem  [11].  Although  there  are  a  number  of  extensions  of  this  theorem 
(see  Theorem  3.1  In  [2]),  we  state  a  simple  version  that  is  sufficient 
for  the  problems  considered  here. 

Theorem  (Trotter-Kato).  Let  s  be  the  generator  of  a  strongly 
continuous  semigroup  e^^  satisfying  |!e'’‘^|i  <  Assume  that  /  is 

d  sequence  of  operators  generating  strongly  continuous  semigroups  e ^ 
satisfying. 

1 )  I  ie'’  !  <  Me’^  ,  N=l,2. . .  . 

11)  -*  AZ  for  z  Z  D.  D  dense  in  z 

i11)  there  exists  Xg  with  Re(  ,\g)  ^  3  such  that  (A-'gl):!  Is  dense 
in  Z.  Then  e  7  -*  e  ^  for  all  z  c  Z  and  the  convergence  is  uniform  for 
t  ^  [O.T]. 

In  terms  of  numerical  analysis,  condition  i)  is  the  stabi 1 i ty  require¬ 
ment  and  condition  ii)  is  the  consistency  requirement.  Thus  convergence 
of  the  approximation  scheme  Is  dependent  upon  having  a  consistent  and 
stable  numerical  scheme  that  also  satisfies  the  technical  condition 
ill).  Any  numerical  scheme  that  does  not  satisfy  conditions  i)  -  iii) 
will  not  produce  an  approximating  sequence  for  the  control  problem  (3.3; 

-  (3.4).  Moreover,  even  if  the  approximating  sequence  is  such  that  i'  - 
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i1i)  are  satisfied,  there  is  no  assurance  that  It  will  be  a  strong 
approximating  sequence  unless  the  operators  and  a'  also  satisfy  i] 

-  iii). 

The  following  convergence  results  are  the  basis  for  numerical 
algorithms  for  computing  sub-optimal  controls.  The  fundamental  results 
are  due  to  Gibson  [3],  Powers  [8]  modified  Gibson's  results  to  apply  to 
Chandrasekhar  algorithms. 


Theorem  3.1.  Let  b^,  be  an  approximating 

sequence  for  the  control  problem  (3.3)  -  (3.4).  Let  n^(t) 
De  the  solution  to  the  Kiccati  equation 

.  -  [/]*  'i(t)  -  -  o" 


n"(t), 


0  >  t 


(3.13) 


wi  th 


n^(T)  »  (3.14) 

i)  If  n(t)  is  the  solution  to  the  Riccati  integral  equation  (3.6), 
then  for  each  t  £  [0,T],  z  s  Z,  r^(t)z  n  (t)2. 

ii;  If  the  approximating  sequence  is  a  strong  approximating 
sequence,  then  n^(t)z  -n(t)z. 
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i1i)  if,  In  addition,  the  control  space  U  is  finite  dimensional 
and  -  vz  for  z  e  Z,  then  for  0  <  t  <  T 


l!^(t)  -  ^(t)||  -  0 


(3.15) 


where 

X  {t)  Is  the  gain  operator  defined  by  (3.7)  -  (3.8)  and.‘.'^(t)  is  the 
solution  to  the  Chandrasekhar  equations 


^/(t) . .  [At)]'  .“(t) 

^  t~(t)  .  -.'•(t)  [/  -  ,V!t)l  , 


0  ^  t  <  T 


(3.16) 


(3.i/) 


w1  th 


L  ^( t)  «  and  <^(T)  =  0 


(3.18) 


A  proof  of  this  theorem  will  appear  in  a  forthcoming  paper.  How¬ 
ever,  it  is  clear  that  from  a  practical  point  of  view  it  is  desirable  to 
have  uniform  convergence  of  the  gain  operator  (i.e.  (3.15))  rather  than 
weak  convergence.  Therefore,  considerable  effort  should  be  devoted  to 
constructing  strong  approximating  sequences  that  have  the  additional 
property  that  .  Such  schemes  are  needed  to  ensure  uniform  con- 
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vergence  of  the  approximating  gain  operators  ^ ^(t)  (obtained  form  a 
finite  dimensional  Chandrasekhar  algorithm)  to  the  optimal  gain  operator 
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IV.  NUMERICAL  PROCEDURES 


We  have  seen  that  the  dynamics  of  the  physical  system  can  be 
formulated  as  an  abstract  model.  A  control  problem  has  been  stated  and 
Its  'solution'  given  without  introducing  any  approximations.  In 
particular,  for  the  linear  quadratic  problem  of  Section  III  one  knows 
that  the  optimal  control  can  be  given  in  feedback  form  as; 


M^(t)  •  Kj(t)9(t)  +  K2(t)ai{t) 

+  K3(tJn(t)  +  K^(tK(t) 

(4.1) 

L 

+  /  K5(t,s)  u  (t.s)  ds 

0  5  XX 

+  /'■  Kg(t,s)  U.(t,s)  dS 
0 


The  time-varying  gains  K^,  Including  the  'functional-gains'  and  Kg, 
are  computed  from  the  ‘solutions'  of  the  Chandrasekhar  Equations  (3.7) 
and  (3.8). 

Calculation  of  these  gains  requires  approximation  of/n  and  a*  as 
discussed  In  Section  III.  In  this  section  we  shall  discuss  some  details 
for  constructing  these  approximations.  Since  for  our  problem  a  •  -  a* 
one  need  only  consider  the  task  of  approximating  a.  It  should  be  noted 
that  this  simplification  is  not  generally  possible. 
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Recall  that  our  state-space  Is  Z  -  R**  x  x  r.^  and  it  is  clear 

that  this  space  is  infinite-dimensional  because  of  ther,^  "functions"  in 
the  Zg  and  Zg  coordinates.  One  useful  way  to  think  about  constructing 
the  approximates  of is  to  Imagine  finite-dimensional  subspaces  Z 
in  which  the  functional  coordinates'  are  approximated  as  linear  com¬ 
binations  of  elementary  functions.  The  choice  here  is  to  use  splines 
[10]. 

Specifically,  we  take  the  interval  [0,1]  and  divide  it  into  N  egual 

subintervals  with  (N+1)  'knots'  at  x-0,  1/N,  2/H . 1.  Consider  a  set 

of  Interpolating  splines  with  the  property  that  { j/n;  «  f.j.  The 
degree  of  these  splines  Is  not  yet  specified;  1n  practice  one  might  use 
piecewise  linear  (hat  functions)  or  piecewise  cubic  elements. 

Define  ef  ^  z,  1  •  0,1,. . . ,N,  by 

ef  *  (0,  0,  0,  0,  4)^(x/L),  0) 

i 

and  for  j  =  1,2 . N  let 


(0, 


(0). 


;  "(i: 

r  J 


0, 


(x/L 


/ 


Finally  let  e_^2“  d*  0.  0.  0.  0,  0)  and  e^^  -  (0,  1,  0,  0,  0,  0) 

The  set  -  ;  e_^^  )C  Z  Is  linearly  independent  and 

hence  its  span  is  a  (2N+3)  dimensional  subspace  of  Z  (say  span  »  Z^). 


For  ease  of  exposition  let  denote  the  set  *  {0,1,2,. ..  ,N}  and 
denote  the  set  ( 1,2,3,. .. ,N). 

Note  that  the  basis  has  been  arranged  so  that  each  e^'^  satisfies 
the  boundary  conditions  needed  to  be  in  o(a).  The  ease  of  this  con¬ 
struction  1s  a  favorable  result  of  our  choice  of  state-space.  That  is, 
one  could  have  eliminated  the  coordinates  and  (n  and  -; )  by 
Id  nt1fy1ng  n(t)  *  (J^  (tj.)  and  E,(t)  =  U^^(t,L).  However,  it  would 

then  be  necessary  to  Incorporate  equations  (1.27)  and  (1.28)  as  boundary 
conditions  in  the  domain  of  a.  This  greatly  complicates  the  task  of 
constructing  basis  elements. 

Continuing  with  the  numerical  aspects  we  shall  display  a  matrix 

representation  of  ^  in  terms  of  the  basis  E^.  Let  ^  :  Z-»  Z^  denote 

N  N  N  N 

the  orthogonal  projection  onto  Z  and  (formally)  write  a  a  p  . 

th  N 

The  J  column  1n  the  matrix  representation  of  will  be  the  image  of 

the  basis  element  under  represented  In  terms  of  the  basis  E^. 

N  N 

since  e  »=  ^  -^e,  the  representation  amounts  to  solving  the  usual  normal 
equations  (see  [10])  for  the  best  approximation  of  -^e  in  terms  of  the 
basis  E^. 

Straight  forward  calculations  reveal  that 

^  i  6^,  ■  col  Ll»0»0t0,0,  0]  , 
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0 


(0) 


A  e 


N  , 
1 


i  e  N 


0 


El 


(1) 

(1) 


0 


-El 

oL' 


and 


j  c 


1 


0 

0 

0 

0 

1 

0 


Note  that  we  have  supressed  notation  Indicating  the  dependence  of 

the  spline  functions  on  the  grid  site  parameter  N.  Also  observe  that 

one  could  approximate  the  second  derivative  operator  (e.g.  by  dlffer- 

enclng  the  first  derivative)  and  hence  allows  to  operate  on  functions 

2  2 

not  smooth  enough  to  be  in W  ’  .  This  would  be  needed,  for  example,  if 
the  were  linear  splines. 
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The  normal  equations  are  of  the  form 


(4.2) 


N  N  N  N 

where  6  Is  the  gram  matrix  with  (i,j)  element  e j  ,  e^-  '  and  R  is  the 

N  N 

matrix  with  (i,j)  element  <6^  ,  ;i  e  j  >  .  From  the  form  of  the  basis 
vectors  it  is  readily  seen  that: 

i)  e^2  e^^are  each  orthogonal  to  the  remainder  of  ; 


i1) 

<  <!2  .  e  ^2  ^  ^ 

=  1/2  ; 

iii) 

<  fiN  N  >  . 

Ia  . 

■  ■  1 

2 

i  v) 

each  e^  is  orthogonal 

to 

each  efl^j  ; 

Hence 

is  block  diagonal  and 

,  In 

particular 

1 

6 

- 

0 

*  1/2 

N 

S 

5 

G^ 

Direct  calculations  reveal  that 
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G 


N 

1 


(i  .j) 


<e^ 

i 


-  £IL  b.{s) 
0  ’ 


ds 


1  I  j  £  • 


and  that 


oL"  (0)  (0) 

L  dP 


+  Me  0.  (1)  $.  (1) 

o  *  V 


*  0i'  (1)  (1) 

■*■  Q  ds^ 


i,j  e  Nj 


(4.3) 


Note  that  the  choice  of  B-splines,  which  are  non-zero  over  only  a 
'small'  portion  of  the  unit-interval,  means  that  the  matrices  will  be 
banded.  The  band-width  increases  with  the  degree  of  the  splines. 


43 


N 

Assembly  of  the  R  matrix  proceeds  in  much  the  same  way  with 
>A  •  It  Is  easily  shown  that 


R  «  1/2 


00 


»! 


(4.4) 


wi  th 


To  l" 

L° 


(  j  e  N.) 


(i.J)  ’  £![--;  (1)  -fy'd)  -  (0)  .y(0) 


’  *1 '(s)  ty  (s)  ds  ] 

0 


i  I  J  £  ^  ^ 


and  =  LR^]  . 

2  1 

N  N 

The  block  structure  of  G  and  R  suggest  a  compatible  partitioning 
N 

of  A  in  the  normal  equations.  It  is  readily  seen  that 
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Formally  then  for  any  N  the  system  {,2.1)  Is  approximated  by 


X  (t)  »  A  X  (t)  +  5  M^(t) 


where  is  as  above  and  5^  is  found  by  projecting  b  onto  Z^. 
The  approximating  solution  is  found  from  the  coordinates  as 


z  {t)  =  I  x.{t)  e^ 


N  N 

The  vectors  x  and  z  are  each  of  dimension  (2N+3).  The  calculations 

N  N 

necessary  to  generate  a  are  easily  done  on  a  computer  once  the 

functions  are  selected. 

In  order  to  illustrate  these  ideas  a  FORTRAN  code  was  constructed. 
The  code  assembles  the  matrices  required  to  evaluate  the  Chandrasekar 
equations  (3.16)  -  (3.18).  Since  many  of  the  matrices  can  be  quite 
sparce,  the  code  employs  partitioned  versions  so  that  most  of  the  zero 
entries  are  not  involved.  For  illustrative  purposes  the  basis  functions 
are  'linear  splines'  [IG]  and  all  physical  constants  are  taken  to  be 
unity.  The  output  operator  '/  in  equation  (3.2)  "reads  out"  the  first 
four  state  components  (i.e.  i{t),U^),  n(t),  n(t));  it  "samples" 

Zg(  t)  =  u^^(t,x)  at  X  =  .2  and  x  =  .4;  and,  it  "samples"  Zg(t)  = 
(u^(t,x)  X  u.(t))  at  X  =  .3  and  x  =  .7.  Thus,  has  rank  8. 
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The  differential  equations  (3.16)  and  (3.17)  were  numerically 
Integrated  backwards  from  the  initial  condition  (3.18)  for  10  time 
units.  It  was  observed  that  the  dependent  variables  k  and  l  were 
changing  very  slowly  at  that  point  and  so  these  were  taken  to  be 
"steady-state"  values.  Shown  in  Figure  3  are  graphs  of  the  kernel 
function  K5(x)  [cf  equation  3.9]  for  N  =  4  and  N  =  8.  These  results  are 
somewhat  preliminary.  In  practice  one  would  continue  to  increase  N  at 
least  until  the  observed  differences  were  "small".  Our  research  code 
had  modest  dimensions  fur  the  various  arrays  and  so  the  maximum 
permissible  N  was  rather  limited.  Note  that  with  N  »  8  the  matrices 
and  in  equations  (3.16)  and  (3.17)  are  1  x  19  and  8  x  19, 
respectively. 
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V.  THE  TWO  BEAM  STRUCTURE 

In  this  section  we  develop  the  equations  of  motion  of  the  structure 
shown  in  Figure  4.  It  is  assumed  that  members  at  the  top  and  bottom  are 
rigid  bodies  connected  by  two  flexible  beams  of  identical  length  L,  the 
structure  is  allowed  to  pivot  about  a  fixed  pivot  point,  and  that  the 
motion  is  In  a  plane.  The  development  of  the  model  is  achieved  by  the 
same  method  used  to  derive  the  equations  of  motion  for  the  single  beam 
with  a  tip  mass.  Since  the  derivation  is  very  similar  to  the  beam-tip 
mass  problem,  we  shall  simply  summarize  the  equations  below. 

Let  u  (t,x  }  and  u  (t,x  )  denote  the  position  of  the  mass  element 
dm  with  respect  to  its  undeflected  position  on  beam  1  and  beam  2, 
respectively.  The  angle  i^;(t)  denotes  the  rotation  of  the  structure 
about  the  pivot  point  and  i|>(t)  is  the  angle  of  rotation  of  the  top  body 
measured  from  its  undeflected  position.  Let  6  denote  the  fixed  angle 
that  each  beam  makes  with  the  bottom  rigid  body  and  note  that  (see  the 
Free  Body  Diagram  in  Figure  5) 

sin  I?  =  h/L  ,  cos  *?  =  (A-C)/2L.  (5,1) 

Using  the  Free  Body  Diagram  in  Figure  4  we  obtain  the  following  set  of 
equations. 

f  sin  e  +  t  =0  (5.2) 

■'4  ^^x 

f^  cos  5  +  fp  =0  (5.3) 


Bottom  Rigid  Body 


-f  cos  c  +  f  sin  e  +  f  cos  e 

xi  yi  X4 


-f  sin  e  -  f  cos  e  -  f  sin  e  - 
^  ■''1  ^’4 
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(5.4) 


*  A/2  (l'yj-fy^)cos  6  +  M^  + 

=  IacmV  (t) 


earn  1 


f 

X 


1 


+  f 

X2 


-m.  Ah  ,.(t) 


(5.5) 


(5.6) 


"2  -  ‘'l  ♦  ■  till  -  f  (A-OjDj)  V  (t)  -  I.jKV  (t) 

^  /SxV„(t.xl)dxl  (S.7, 


Beam  2 


f  *.  f  .  .,  Ah  . 
-3  'X4  2  B2  '*' 


(5.8) 


^3  ^  \  “  [(f  -  7  (^)  V  (t)  -  ^-(t)  +  u2^^(t,^-)j 


9) 


^'3  -  M4  ■  *  CI24  -  |(A-C)n,2]  V  (t)  -  l24Kr  (t) 


+  A  ,Jx^ij^^^(t,x^)dx^ 


(5.10) 
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Top  Rigid  Body  (mass  »  m^) 

(f  -f  )cos  0  +  (f  -f  )s1n  0  =  m  h(Ki’(t)-V  (t))  (5.11) 

3  ^2  ^2  ^3  ^ 


0  +  ify^*fy^)cos  0  -  0  (5.12) 

=  ICCM  V  (t)]  (5.13) 

Observe  that  equations  (5.4),  (5.11)  and  (5.13)  can  be  combined  to 
produce  the  algebraic  equation 


/COS0 

'm^F' 


CK 

^^CCM 


) 


>3  "2 


A  sin  e 
ACM 


(f 


X 


1 


/ 


=  ( 


sine 


CK  V 

71 — ' 

‘^^CCM 


^3  ^2 


+ 


A  cos  0 
ACM 


(5.14) 


ni!<i 

^ACM 


(Ml  +  %) 


TCM 


(M-  +  M3)  - 

‘acm 


'AC 


Let  denote  the  column  vectors  =  col  (f^  ,f^  ,f^  .f^  )  = 


respectively.  Combine  the  algebraic  equations  (5.5)  (5.8)  (5.12)  and  (5.14) 
into  the  one  linear  equation 


so 


We  note  that  M  exists  and  therefore  7^  can  be  determined  as 
a  function  of  f^  and  R.  In  order  to  complete  the  model  we  need  to 
provide  equations  for  the  elastic  motion  of  the  beams  and  the  corresponding 
boundary  conditions.  These  equations  are; 


.(t.xM 


•  x^ri;‘(t)  -  kV  (t)] 


(5.19) 


u^^^(t,x^) 


.1 


^  (t)  -  kV  (t)] 


(5.20) 


with  boundary  conditions: 

at  x^  •  0  f^^  » 

u^(t,0)  »  0 
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!M2  = 

f  m  E^I^u^xxx^^’**^  (5.22) 

^2 

u^(t,L)  -  0 

u^^(t.L)  »  (l-K)A(t) 


"4  ■ 

u^{t.O)  -  0 

u^y(t,0)  =  K(}»(t) 

M3  •  E2l2,2^^(t,L) 

4t  x'  L  \ 

/  U^(t,L)  «  0 
\  u^j^(t,L)  »  (1+K)0(t) 


(5.23) 


(5.24) 


Let  z(t)  denote  the  vector 

2(t)  =  («(t).^*(t).  U^^^(t,.).0(t)J(t).  u2^^(ti).  u2,.(t..)^ 

and  note  that  z(t)  [R^  x  (0,L.  Therefore,  the  state  of 

the  system  1s  described  by  four  functions  d),  d),  \i),  ij;  and  four  L2 

valued  function  ij^  ,  u^^,  u^  ,  u^ 
xx’  “  t’  xx’  t 
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Equations  (5.19)  (5.20),  (5.4)  and  (5.11)  can  be  written  as  the 
first  order  system 

Fi(t)  •  H2(t)  +  G  (5.25) 

2  2 

where  H  is  a  boundary-differential  operator  on  2  =  [R  x  L2  (0,L,  F  )] 
2  2 

X  [R  X  L2  (0,L,  F  ],  and  F  is  a  bounded  non-singular  operator.  Let 
A  »  (F"^H)  8  «  F"^  G  (5.26) 

and  define 


0 

(A)  =  (z  c  Z 

/  2  satisfies 

Hi) 

(5.27) 

where  z  »  (z 

1’  ^2’ 

^6 

.  Zj(.),  2g(.)) 

and  Hj)  i) 

Z3(.).  Z4(. 

).  Zg 

,(.) 

all  belong  to 

’^(0,L). 

ii) 

0 

n 

0 

♦SI 

(0) 

iii) 

Z4(U  =  0 

’  ^8 

(0) 

=  K  Zg 

iv) 

Zg(0)  =  0 

•  ^8 

(0) 

v) 

Zg(L)  -  0 

’  ^8 

(L) 

=  (1+K)  Zg. 

(5.26) 

The  specific  form  of  F,  H,  G,  A  and  B  can  easily  be  constructed 
from  equations  (5.4),  (5.11),  (5.19)  and  (5.20)  by  combining  equation 
(5.15)  with  the  equations  (5.6),  (5.7),  (5.9),  and  (5.10).  As  in 
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the  problem  with  the  single  beam,  it  is  necessary  to  show  that  the 
operator  a  defined  by  (5.26)  -  (5.28)  generates  a  dynamical  system 
on  Z  .  Numerical  algorithms  will  be  similar  to  the  one  In  Section 
4  above. 


55 


REFERENCES 


1.  Adams,  R.  A.,  Sobolev  Spaces.  Academic  Press,  New  York,  1975. 

2.  Banks,  H.J.,  Burns,  J.A.  and  Cliff,  E.M.,  Parameter  estimation  and 
Identification  for  systems  with  delays,  SIAM  ih  Control  and 
Optimization,  1981,  791-828. 

3.  Gibson,  J.S.,  "The  Riccati  integral  equations  for  optimal  control 
problems  on  Hilbert  spaces,"  SIAM  J.  Control  and  Optimization,  17, 
1979,  533-565. 

4.  Hllle,  £.,  and  Phillips,  R.S.  Functional  Analysis  and  Semi-groups, 
AMS  Colloq.  Pub.  Vol.  31,  1957. 

5.  Meirovitch,  L.  Elements  of  Vibration  Analysis,  McGraw-Hill,  New 
York,  1975. 

6.  Miele,  A.,  Flight  Mechanics,  Addison  Wesley,  Reading,  MA,  1962. 
pp.  28-34. 

7.  Pazy,  A.  Semigroups  of  Linear  Operators  i'nd  Appl  1  cat  ions  to 
Partial  Differential  Equations .  Springer-V.'rlag," ^ew  Vor^ ,T983. 

8.  Powers,  R.K.,  "Chandrasekhar  equations  for  distributed  parameter 
systems,"  Ph.D.  Dissertation,  VPI&SU,  Blacksburg,  VA,  1984. 

9.  Rudin,  W.  F.,  Principles  of  Mathematical  Analysis .  McGraw-Hill,  New 
York,  1976. 

10.  Schultz,  M.H.,  Spline  Analysis,  Prentice-Hall,  Englewood  Cliffs, 

NJ,  1973. 

11.  Trotter,  H.F.,  "Approximations  of  Semi-groups  of  operators," 

Pacific  ^  Math.  8,  1958,  887-919. 

12.  Walker,  J.A.,  Dynamical  Systems  and  Evolution  Equations ,  Plenum 
Press,  New  York,  1980. 


56 


X 


Mass  A 
{main  frame) 
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Figure  3:  Functional  Gain  Kd») 


Pointing  Axis 


PART  II 


EFFECTS  OF  TIME  DELAY.  ACTUATOR  DYNAMICS 
AND  SYSTEM  DAMPING  CN  THE  FEEDBACK 
CONTROL  OF  A  FLEXIBLE  CABLE 
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An  Example  of  the  Effect  of  Time  Delay  on 
Systems  with  Feedback  which  are  Governed  by 
Partial  Differential  Equations 


I.  Introduction 

The  purpose  of  this  work  Is  to  offer  some  insight  into  some 
problems  associated  with  the  design  of  controllers  for  systems  governed 
by  partial  differential  equations.  The  method  of  analysis  is  based  on 
the  concept  of  transfer  functions  which  relate  the  output  or  response  of 
a  system  to  the  input.  For  simple  structures,  the  displacement 
(velocity,  or  slope)  can  be  related  to  the  force  applied  by  an  exact 
transfer  function.  The  accuracy  of  the  transfer  function  Is  limited 
only  by  the  accuracy  of  the  mathematical  model  used  to  describe  the 
structure.  This  open-loop  transfer  function  can  be  used  to  construct 
the  closed-loop  transfer  function 'for  the  case  where  the  displacement 
(velocity,  or  slope)  Is  fed  back  to  the  force.  By  examining  the 
characteristic  equation  of  the  closed-loop  transfer  function  the 
stability  of  the  closed-loop  system  can  be  analyzed.  This  method  will 
be  applied  to  a  flexible  cable  under  tension. 

It  should  be  pointed  out  that  although  only  a  simple  cable  mass 
element  will  be  considered  here,  the  procedure  can  be  used  to  examine 
more  complex  structures  [l].  The  advantage  of  the  formulation  used  here 
is  that  it  provides  information  In  the  form  of  root  locus  plots  which 
are  familiar  to  most  and  hence  adds  to  the  physical  interpretation  of 
the  results.  The  inclusion  of  time  delay  in  the  problem  does  not  change 
the  method  of  analysis. 
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It  Is  likely  that  in  the  near  future  almost  all  control  systems 
will  be  digital  in  nature.  Consequently  the  control  command  will  always 
respond  to  the  sensor  inputs  with  a  delay  of  at  least  one  interval  of 
sampling  time  associated  with  the  controller.  The  contention  to  be 
demonstrated  here  is  that  regardless  of  how  small  this  delay  may  be,  it 
can  cause  some  instability  in  systems  with  feedback.  Because  of  its 
simplicity  and  familiarity,  the  analysis  will  be  applied  to  the  cable 
problem  as  indicated  previously. 


II.  Flexible  Cable  Analysis 

The  vibrating  cable  is  considered  to  be  a  continuous  or  distributed 
parameter  system,  that  is,  one  governed  by  a  partial  differential 
equation.  As  shown  in  Figure  la,  f(x,t),  P(x),  and  T(x)  are  the 
distributed  force,  mass  density,  and  tension  in  the  cable,  respectively, 
expressed  as  a  function  of  position,  x  along  the  cable.  For  this 
analysis,  negligible  structural  damping  and  no  transverse  stiffness  is 
assumed. 

The  equation  of  motion  describing  the  transverse  motion  of  the 
cable  can  be  obtained  by  examining  a  differential  element  of  the  cable. 
Figure  lb  represents  the  free  body  diagram  corresponding  to  a 
differential  element  of  cable  of  length  dx.  Applying  Newton's  second 
law  in  the  vertical  direction,  assuming  small  deflections,  only  vertical 


64 


motion  and  Ignoring  second-order  terms  In  dx,  we  find  that  the  governing 
partial  differential  equation  of  motion  of  the  cable  Is  given  by  [2] 

^  CT(x)  +  f(x,t)  -  o(x)  .  0  <  X  <  L  (2.1) 

3t 

While  Eq.  (2.1)  Is  the  general  equation  of  motion  of  the  cable,  It  can 
be  simplified  by  certain  appropriate  assumptions.  In  addition,  any 
solution  of  Eq.  (2.1)  will  depend  on  the  particular  boundary  conditions 
of  the  cable  configuration  under  consideration. 


Fixed-free  cable  with  a  discrete  mass  at  end 

The  configuration  to  be  Investigated  Is  a  cable  of  length  L  with 
constant  mass  per  unit  length,  o,  subjected  to  a  constant  tension  T. 

This  cable  Is  fixed  at  one  end,  while  the  other  end  Is  free  with  a 
concentrated  mass  attached.  A  control  force  Is  acting  vertically  on  the 
mass  at  the  free  end.  Furthermore,  there  Is  no  distributed  force  f(x,t) 
acting  along  the  cable,  (see  Figure  2a).  For  this  particular  problem 
Eq.  (2.1)  reduces  to 

P  ^  j  ,  OlxlL  (2.2) 

3 3x^ 

with  the  associated  boundary  condition  at  the  fixed  end  x  ■  0, 

y(O.t)  -  0  (2.3) 

By  writing  Newton's  second  law  for  the  free  end,  shown  in  Figure  2b,  the 
boundary  condition  for  the  free  end  x  =  L  becomes. 


m  L"y..(x.U 

3r 


x-L 


(2.4) 
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Eq.  (2.2)  can  now  be  solved  for  the  lateral  deflection  y{x,t). 

In  the  course  of  obtaining  this  solution  we  seek  a  relationship  between 
the  velocity  at  the  end  (x-L),  (output)  due  to  a  force  at  the  same  end 
(input).  This  relationship  Is  best  expressed  In  terms  of  a  transfer 
function  which  can  be  obtained  by  taking  the  Laplace  transform  of  eq. 
(2.2)  with  respect  to  the  time  variable  and  evaluating  the  result  at 
x*L.  Since  we  are  interested  In  the  velocity,  the  result  must  be 
multiplied  by  the  Laplace  variable, s  which  Is  equivalent  to  the  time 
derivative.  The  desired  transfer  function  Is  given  by 

y(L.s)  .  _ slnh  /ofT  sL _ 

f|^  >!pT  cos^l^/^  sL  +  ms  s^nh  sL  (2.5) 

In  order  to  simplify  the  problem  further,  we  will  set  p“T»l,  In 

addition  it  is  convenient  to  allow  the  end  mass  to  go  to  zero  since  this 

term  only  effects  the  uncontrolled  frequenclas  and  not  the  system 

stability.  Equation  (2.5)  now  reduces  to 

y(l-.s)  ,  slnh  s  .  G(s)  (2.6) 

cosh  s 

Equation  (2.6)  Is  the  open  loop  transfer  function  relating  the  velocity 
at  the  end  to  the  force  at  the  end  of  the  fixed-free  cable. 

With  eq.  (2.6)  as  a  starting  point  we  can  now  examine  the  system 
characteristics  for  both  open  and  closed  loop  control  using  root  locus 
analysis. 
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■  ^  ii  ~ 


^  ^ 


III.  Flexible  Cable  -  Stability  and  Control 

The  system  natural  (open  loop)  frequencies  can  be  obtained  by 
setting  the  denominator  of  the  transfer  function  equal  to  zero  and 
solving  for  s,  where  s  »  n  +  iw.  It  Is  easily  verified  that  the 
solution  of  cosh  s  ■  0  from  Eq.(2.6)w1ll  yield  Imaginary  values  of  s 
corresponding  to  the  frequencies  associated  with  a  fixed-free  cable. 
These  are  given  by 

n  -  0.  01  =  21^  ^  k  -  ±1,±2,...  (2.7) 

where  k  is  an  Integer  which  serves  to  Identify  a  frequency. 

Feedback  Control 

The  system  has  open  loop  Input,  fj^  and  output  y  related  through 
G(s).  We  now  close  the  loop  by  feeding  back  velocity  with  a  gain  K,  see 
Figure  3. 

The  closed  loop  transfer  function  can  be  shown  to  be 

G(s)^.l  =  [1  +  KG{s)]‘^  G(s)  (3.1) 

The  closed-loop  system  dynamics  can  be  determined  from  the  closed-loop 
characteristic  equation 

1  +  K  G(s)  =  0 

I  Now  all  the  classical  tools  for  control  design  can  be  used  [3]. 
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(3.2) 


For  velocity  feedback  at  the  cable  end,  the  open-loop  transfer  function 
was  shown  to  be 


G(s) 


sinh  s 
cosh  s 


(2.6) 


Using  Eq.  (3.2)  the  associated  closed-loop  characteristic  equation 
becomes 


cosh  s  +  K  sinh  s  =  0  (3.3) 

If  we  let  s  =  n  -t-  iw  we  can  obtain  an  analytic  solution  for  the  feedback 
by  velocity.  The  results  are  given  by 

Tanh  n  *  -K  j  w  *  (m  odd) 

(3.4) 

Tanh  n  =  ^  ^  (m  even) 


The  root  locus  plot  is  shown  in  Figure  4. 

It  is  seen  from  the  root  locus  that  all  the  modes  are  controlled  by 
feeding  back  the  velocity  at  the  end  of  the  cable  to  a  force  at  the  end. 
Furthermore  no  special  filter  is  needed  to  process  the  sensor  signal. 

For  this  particular  case  all  modes  are  affected  the  same  for  a  given 
gain.  For  gains  less  than  one,  the  frequencies  of  vibration  are  the 
same  as  the  open-loop  frequencies  but  the  motion  is  damped  out.  For 
gains  greater  than  one  the  frequencies  jump  to  those  of  a  fixed-fixed 


68 


cable  and  the  motion  is  also  damped.  At  a  gain  of  one,  n  ■  -  «  and  1t 
can  be  shown  that  the  system  comes  to  rest  in  a  finite  time  [4]. 

Although  this  simple  problem  can  be  solved  analytically,  the  more 
general  problem  to  follow  would  become  analytically  tedious. 
Consequently,  a  numerical  method  is  used.  Host  computer  routines  that 
generate  root-loci  need  to  be  supplied  with  polynominal  functions  or 
factors  [5].  In  this  work  the  solution  of  a  transcendental  equation  is 
required.  One  method  of  solution,  for  the  generation  of  the  root  loci, 
is  the  use  of  a  computer  code  which  includes  an  IMSL  routine  called 
ZSCNT  which  solves  for  the  roots  of  a  set  of  non-linear  simultaneous 
equations  [6].  The  complex  characteristic  equation,  Eq.  (3.3)  is  one  of 
these  sets  when  separated  into  real  and  imaginary  equations.  The 
procedure  is  to  apply  repeatedly  the  ZSCNT  routine  starting  with  the 
open-loop  pole  position  with  zero  gain.  Then,  incrementing  the  gain, 
solve  for  the  first  closed-loop  pole  position.  This  pole  is  then  used 
as  the  start  point  for  solving  for  the  pole  at  the  next  higher  gain 
value  and  so  on.  This  method  is  found  to  be  highly  reliable  for 
calculating  root-loci,  although  occasionally  sensitive  to  abrupt  changes 
in  locus  direction.  A  program  listing  with  an  example  input  and  output 
file  is  given  at  the  end  of  this  report. 


Feedback  Control  with  Time  Delay 

Time  delay  is  now  introduced  into  the  system  such  that  the  feedback 
of  velocity  is  described  by  the  d1ag''am  in  Figure  5.  As  we  have  deter- 
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mined  previously,  the  open-loop  velocity  force  transfer  function  for  the 
fixed-free  cable  Is 


G(s)  - 


sinh  s 
cosh  s 


(2.6) 


Adding  the  time  delay  of  e"^^  to  the  feedback  loop  yields  a  closed  loop 
transfer  function  of  the  form 


G(s)cl 


G(s) 

1  +  Ke*^^  G(s) 


which  results  In  the  characteristic  equation 


(3.5) 


cosh  s  +  K  e’^ '  sinh  s  »  0 


(3.6) 


or 


cosh  s  +  K  (cosh  ST  -  sinh  sx  )  sinh  s  =  0  (3.7) 

By  letting  s  =  n  +  1^,  and  Incrementing  the  gain,  K  for  a  fixed  time 
delay,  t  ,  the  root  locus  for  the  cable  with  delayed  feedback  of 
velocity  can  be  obtained.  The  root  locus  for  time  delays  of  0.1,  0.5, 
and  0.9  are  shown  in  Figures  6,7,  and  8,  respectively.  These  diagrams 
show  clearly  that  unlike  the  no  delay  case  that  for  any  delay  there  are 
associated  unstable  roots.  A  relationship  between  the  delay  time  and 
the  unstable  root.,  can  be  determined  as  follows: 
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Starting  with  the  characteristic  equation  for  the  cable  with  delay 
eq.(3.7)  we  can  note  that  it  has  the  form  of 

F  (s,  k,  T  )  •  0  (3.8) 

Letting  s  =  n  +  ioj  and  expanding  yields  a  real  equation  and  an  imaginary 
equation  of  the  forms 

R  (n,  <1.,  K.  T  )  =  0  (3.9) 

and 

I  (n,  ij,  K,  T  )  »  0  (3.10) 

If  we  select  a  particular  value  of  t  we  can  suppress  its  dependence  in 
eqs.  (3.9) and  (3.10).  Further  we  can  note  that  the  real  and  imaginary 
parts  of  the. solution  along  any  branch  of  the  root  locus  are  functions 
of  the  gain  K.  Hence  eqs.  (3.9)  and  (3.10)  take  the  form 

R  (n(K),  4K),  K)  =  0  (3.11) 

and 

I  (n(k),  4K),  K)  =  0  (3.12) 
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By  using  the  derivative  chain  rule  we  find 
3R  dn  ,  d(i)  ,  ^  « 

3n  dK  3w  HI?  3K  ”  (3.13) 


and 


31^  dn 

3n  HR  3oj  HI?  3l?  ■ 


(3.14) 


We  are  Interested  in  solving  for  the  quantity  ^  at  the  point  where 

dK 

K  =  n  a  0.  Evaluated  at  these  conditions  describes  how  the  real 

OK 

part  of  the  solution  departs  from  the  zero  gain  roots  on  the  imaginary 
axis.  A  positive  value  would  yield  instability  for  a  small  gain. 
Solving  for  4!!  or  n'  yields 


n  = 


R  I 

iij  N  N  'jk) 

R„  I  -  R  L 

n  'ii  u,'  n 


(3.15) 


where  the  subscripts  indicate  the  partial  derivatives  in  (3.13)  and 
(3.14). 

Applying  eq.  (3.15)  to  the  characteristic  equation  for  our  system 
(3.7)  yields  the  simple  result 

n'  =  -  cos  T  (3.16) 
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Since  k  must  be  an  integer,  these  values  describe  particular  frequencies 
which  would  be  driven  unstable  for  small  gain.  Also  from  Figures  6,  7, 
and  8  it  can  be  seen  that  the  cycles  are  repeated  as  indicated  by  eq. 
(3.18). 

IV.  Feedback  Control  with  First  Order  Actuator 


The  previous  investigation  showed  the  effect  of  a  control  system 
which  included  pure  delay  on  the  stability  of  a  fixed-free  cable.  The 
purpose  of  this  section  is  to  show  the  effect  a  first  order  actuator 
might  have  on  the  cable  stability  both  with  and  without  a  pure  delay. 


First  Order  Actuator  without  delay 

We  will  initiate  the  study  by  placing  a  first  order  actuator  in  the 
system  as  described  by  Figure  9.  Here  the  time  constant  associated  with 
the  first  order  system  is  designated  by  t.  From  the  figure,  the  closed 
loop  transfer  function  can  be  shown  to  be 


G(s) 


CL 


G(s)  _ 


(4.1) 


where  G(s)  is  given  by  eq.  (2.6).  The  resulting  characteristic  equation 


cosh  s  +  s  cosh  s  +  K  sinh  s 


C 


(4.2) 


By  letting  s  *  n  +  iu,  and  Incrementing  the  gain  K  for  a  fixec  value  of 
T,  the  root  locus  for  the  cable  with  velocity  feedback  through  a  first 
order  actuator  can  be  obtained  using  procedures  described  previously. 

The  resulting  root  locus  for  values  oft  of  0.1  and  0.5  are  shown 
in  Figures  10  and  11  respectively.  These  diagrams  show  that  the  first 
order  actua'^.rr  by  itself  does  not  lead  to  any  unstable  roots  for  any 
gain.  However  the  strong  stability  characteri tics  displayed  by  direct 
feedback  as  shown  in  Figure  4  are  severely  reduced,  especially  at  the 
higher  frequencies.  This  result  Is  easily  explained  by  noting  that  the 
first  order  actuator  has  a  finite  bandwidth  and  hence  does  not  respond 
very  well  to  the  frequencies  outside  of  this  bandwidth.  Typically  the 
bandwith  of  a  first  order  system  is  characterized  by  uig  »  1/t.  For  our 
two  cases,  'jJg  =  10  and  2  for  Figures  10  and  11  respectively. 
Consequently,  frequencies  above  these  values  are  not  affected  as  much  as 
those  below.  However,  as  indicated  previously,  the  control  over  those 
frequencies  below  the  respective  bandwidth  frequencies  is  considerably 
reduced  over  the  case  where  the  actuato**  dynamics  is  ignored.  One  would 
expect  that  to  control  even  low  frequencies,  an  actuator  with  a  large 
bandwidth  is  required. 

First  Order  Actuator  with  Delay 

We  can  now  include  a  pure  delay  in  conjunction  with  a  first  order- 
actuator.  This  situation  might  occur  if  a  digital  control  system  is 


used  to  drive  the  actuator.  The  block  diagram  displaying  this  type  of 
set  up  is  shown  in  Figure  12.  From  this  diagram  we  can  develop  the 
closed  loop  transfer  function  which  is  given  by 
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(4.3) 


wnich  leads  to  the  characteristic  equation 


cosh  s  +  ST  cosh  s  +  K  (cosh  st  -  sinh  ST)sinh  s  =  0 


(4.4) 


Proceeding  in  the  same  manner  as  before,  the  root  locus  for  a  cable  with 
velocity  feedback  thru  a  first  order  actuator  with  delay  can  be 
determined. 

The  root  loci  for  the  cases  where  t  ••  0.1  and  0.5  are  shown  in 
Figures  13  and  14  respectively.  It  can  be  observed  that  in  both  cases, 
the  pattern  of  stable  and  unstable  roots  is  the  same  as  that  which 
occurred  with  pure  delay.  The  first  frequency  at  which  instability 
occurs  however  is  slightly  lower  for  the  case  of  a  first  order  actuator 
with  delay  than  that  with  pure  delay. 


V.  Feedback  Control  with  System  Damping 


The  various  types  ot  damping  that  could  be  considered  are 
structural,  viscous,  and  Coulomb  [8].  The  viscous  damping  model  is 


r6 


chosen  here,  as  this  model  is  the  most  easily  understood  and  most  widely 
used. 

The  Fixed-Free  cable  model  for  a  constant  tension  cable  with 
viscous  damping  is  given  by  eq.  (2.1)  or 

cy^(x,t)  +oy^^(x,t)  = 

where  the  viscous  damping  is 

cy^(x,t)  »  -f(x,t)  :) 

The  boundary  condition  at  the  fixed  end,  x»0  is 

y(0,t)  -  0  (2,3) 

and  at  the  free  end,  x-L  is 

j  x=L  "  (2.4) 


The  transfer  function  of  this  model  can  be  obtained  by  taking  the 
Laplace  transform  of  eq.  (5.1)  with  respect  to  the  time  variable  and 
evaluating  the  result  at  x»L.  Again  we  are  interested  in  the  velocity 
force  transfer  function.  For  the  case  where  the  end  mass  is  zero,  this 
transfer  function  now  has  the  .'orm 


'^7 


(5.3) 


G(s) 


s  sinh  /s*  +  cs 
+  cs  cosh  +  cs 


The  system  diagram  for  the  damped  system  is  illustrated  in  Figure  3  with 
G(s)  now  given  by  eq.  (5.3). 

The  closed-loop  characteristic  equation  is  given  by 

+  cs  cosh  Vs^  +  cs  +  K  s  sinh  Vs^  +  cs  (5.4) 

By  letting  s  *  n+ioj  for  fixed  damping  coefficient,  c  and  incrementing 
the  gain,  K,  the  root  locus  can  be  obtained  as  before. 

The  addition  of  damping  to  the  closed-loop  system  merely  shifts  the 
root  loci  from  poles  with  zero  real  parts  to  poles  that  have  negative 
real  parts.  The  diagram  will  look  like  Figure  4  with  starting  and 
ending  points  shifted  into  the  left  half  plane  and  again  no  unstable 
roots . 

Feedback  Control  with  System  Dampinq  and  Delay 

A  more  interesting  result  is  obtained  by  adding  time  delay  to  a 
system  with  damping.  The  closed-loop  characteristic  equation  of  this 
case  is  given  by 

rs^  I-  cs  cosh  '  S''  +”  «  +  Ks  (cosh  ST  -  sinh  st)  sinh  -  s*  +  cs  (5.5) 
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By  letting  s  »  n  +  I'^'and  Incrementing  the  gain,  K  for  a  fixed  time 
delay,  and  a  fixed  damping  coefficient,  c,  the  root  locus  ^or  the  cable 
with  damping  and  delayed  feedback  of  velocity  can  be  obtained.  The  root 
locus  for  time  delay  of  0.1  and  damping  coefficients  of  0.1  and  1.0  are 
shown  in  Figures  15  and  16  respectively.  These  diagrams  are  comparable 
to  the  diagram  in  Figure  6,  but  shifted  to  the  left  by  an  amount 
proportional  to  the  damping  coefficient.  This  observation  leads  to  the 
conclusion  that  given  sufficient  damping  the  system  is  stable  for  small 
values  of  gains  even  with  time  delay. 

VI.  Closure 

The  cable  problem  discussed  in  this  report  can  be  made  extremely 
stable  by  feeding  back  velocity  at  the  end  of  the  cable  to  a  cc-located 
force  at  the  end  of  the  cable.  However  it  was  shown  that  if  any  delay 
at  all  is  used  in  implementing  this  feedback  control,  the  system  becomes 
unstable.  This  result  was  determined  analytically  with  no  restrictive 
assumptions  other  than  those  maoe  for  the  original  governing  equation. 
The  results  were  confirmed  for  selected  numerical  calculations. 

In  addition  the  results  of  a  first  order  actuator  and  a  first  order 
actuator  with  delay  were  presentn-*  indicafittg  that  any  delay  leads  to 
instability.  Without  de'.ay  the  first  order  actuator  does  not  lead  to 
ins. ability  but  is  considerably  less  effective  in  stabilizing  the  cable 
than  if  no  actuator  dynamics  were  present.  Further,  the  effec.a  at 
higher  frequencies  are  limited  due  to  the  bandwidth  or  the  actuator. 
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The  addition  of  damping  shifts  the  whole  root  locus  to  the  left  and 


allows  a  system  with  pure  delay  to  remain  stable  for  a  small  amount  of 
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Program  Listing 

C  THIS  ROUTINE  SOLVES  A  SYSTEM  OF  TRANSENDENTAL  NONLINEAR  EQUATIONS 

C 

C  IT  USES  IMSL  ROUTINE  ZSCNT 

C  THE  VARIABLES  USED  ARE:  PAR(l)  »  G  »  GAIN,  PAR(2;  *  TIME  DELAY 

C  X(l)  =  R  =»  REAL  PART  OF  ROOT,  X(2)  =  W  =  IMAGINARY  PART  OF  ROOT 

C  S  *  COMPLEX  {R,W),  K  =  NO.  OF  ROOTS,  NR  =  NO.  OF  POINTS  PER  ROOT 

C  THE  PROGRAM  READS  FROM  DATA  FILE  3  AND  WRITES  TO  FILE  4  . 

C  CHAR  =  CHARACTERISTIC  EQUATION  (THIS  CASE  FIX-FREE  CABLE  W/DELAY 

COMPLEX  SH,2,CH,Sr ,CO,S,CHAR,TS 
EXTERNAL  FCN 

DIMENSION  X(50) ,F(50) , PAR (50, ,WK(50) 

C  READ  IN  TIME  DELAY 

READ  (3,*)  PAR(2) 

C  READ  ONE  AT  A  TIME  EACH  OF  19  ROOTS,  INITIALIZE  GAIN  =  0 

DO  30 

PAR( 1 ) =0 . 0 

READ  (3,*)  X(1),X{2) 

WRITE  (4,5)  X( 1) ,X( 2) ,PAR(2) 

5  FORMAT  ( i0X,33HSTART  REAL  IMAGIN,  AND  TIME  DELAY ,/ , 3El 5 . 4 , // ) 

WRITE  (4,2) 

2  F0RMAT(6X,7HREAL  S,9X,7HIMAG.  S,7X,10HNORM  ERROR , 5X , 4HGAIN ) 

C  8  SIGNIFICANT  FIGURES,  TWO  EQUATIONS  (I.E.  REAL  AND  IMAGINARY) 

NSIG  =>  8 
N  =  2 

UMAX  =  200 
KNN  =  N 

C  SOLVE  EQUATIONS  FOR  100  POINTS  FOR  EACH  ROOT 

DC  20  NR  =  1,100 

CALL  ZSCNT  ( FCN , NS IG , N , ITMAX , PAR , X , FNORM , WK , lER ) 

WRITE(4,1)  (X(I),  1=1, KNN),  FNORM,  FAR(l) 

1  FORMAT  (4E15.4) 

C  INCREMENT  THE  GAIN  SLOWLY  SO  THE  SOLVER  TO  FOI-LOW  THE  ROOT  LOCUS 

20  ?AR(1)  =  1 .04*(PAR( l)+0.02) 

WRITE',  1,4) 

4  FORMAT  . IHi) 


9? 


30 


CONTINUE 
STOP 
END 

SUBROUTINE  FCN  (X,F,N,PAil) 

COMPLEX  SH,Z ,CH,SI ,CO,S,CHAR,TS 
DIMENSION  X(N) ,F(N) ,PAR(1) 

C  ABREVIATE  SIN,  COS,  SINK,  AND  COSH 
SH ( Z ) » ( CEXP ( Z ) -CEXP (-Z))/(2. 0,0.0) 

CH ( Z ) * ( CEXP ( Z ) +CEXP (-Z))/(2. 0,0.0) 
SI(Z)aCSIN(Z) 

CO(Z)«CCOS(Z) 

R-X(l) 

W=X(2) 

G»PAR(1) 

S*CMPLX  (R,W) 

TS  »  S  *  PAR(2) 

C  THE  CHARACTERISTIC  EQUATION  GOES  HERE 
CHAR-CH ( S ) +G* ( CH ( TS ) -SH ( TS ) ) *SH ( S ) 

F( 1)«REAL(CHAR) 

F( 2)»AIMAG{CHAR) 

RETURN 

END 


Data  Input  File 


0 . 1 

0.0  1.570796 
0.0  4.712389 
0.0  7.853982 
0.0  10.995574 


output  For  One  Root 


START  REAL  IMAGIN.  AND  TIME  DELAY 


0.0 

0.1571E+01 

REAL  S 

IMAG.  S 

0.0 

0.1571E+01 

-0.2059E-01 

0.1574E+01 

-0.4211E-01 

0.1578E+01 

-0.6462E-01 

0.1581Et01 

-0.8820E-01 

0.1585E+01 

-0 . 1129E+00 

0.1589E+01 

-0.1389E+00 

0,1593E+01 

-0 . 1663E+00 

0.1598E+01 

-0 . 1952E+00 

0.1603E+01 

-0 . 2257E  +  00 

0.1609E+01 

-0 . 2580E+00 

0.1615E+01 

-0 . 2924E+00 

0.1621E+01 

-0 . 3291E+00 

0.1629E+01 

-0 . 3685E+00 

0.1638E+01 

O.lOOOE+00 


NORM  ERROR 

GAIN 

0.9858E-13 

0.0 

0 .5090E-12 

0.2080E-01 

0.1089E-13 

0 . 4243E-01 

0 .2269E-11 

0 . 6493E-01 

0.9173E-12 

0 . 8833E-01 

0 . 2791E-12 

0 . 1127E+00 

0.9437E-13 

0. 1380E+00 

0 . 1360E-11 

0 .1643E+00 

0.2402E-11 

0.1917E+00 

0 .1796E-11 

0.2201E+00 

0 .9998E-13 

0.2497E+00 

0 . 1628E-11 

0 . 2805E+00 

0.1374E-11 

0.3125E+00 

0 . 3695E-12 

0 . 3458E+00 
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